#include "Dwarf.Unmanaged.Omp.h"

Solver :: Solver()
{
    length = 0;
	Wlength = -1;
	W = NULL;
}

//Dispose arrays
Solver :: ~Solver()
{   
    for(int i = 0;i < length; i++){
        delete[]complexArray[i];
    }
    delete []complexArray;
	delete W;
}

void Solver :: init()
{
    this->complexArray = new complex<double>*[length];
    for(int i = 0; i < length; i++){
        this->complexArray[i] = new complex<double>[length];
    }
}

//Transposition of matrix
void Solver :: transpose(complex<double>** matrix)
{
    complex<double> temp;
    for (int i = 0; i < length; i++)
    {
        for (int j = i + 1; j < length; j++)
        {
            temp = matrix[i][j];
            matrix[i][j] = matrix[j][i];
            matrix[j][i] = temp;
        }
    }
}

//1-dimension Fast Fourier Algorithm
void Solver::FFT(complex<double> *complexLine) 
{
    //variables for the FFT
    int bitwidth = (int)(log((double)length)/log(2.0));
    complex<double> w,t,u,v;
    //the complex array splits to real array and complex array 
	if (W == NULL || Wlength != length) 
	{
		Wlength = length;
		W = new complex<double>[length];

		//compute phase weights
		for (int i = 0; i < length; i++)
		{
			W[i].real(cos(-((2* i * M_PI) / length)));
			W[i].imag(sin(-((2* i * M_PI) / length)));
		}
	}

    //reorder for butterflies
    int firstBit = length >> 1;
    int index = 0;
    int reverseIndex = 0;
    for (;;) {
        if (reverseIndex < index) {
            w = complexLine[index];
            complexLine[index] = complexLine[reverseIndex];
            complexLine[reverseIndex] = w;
        }
        index++;
        if (index == length)
            // Break here, because the code below doesn't terminate for all 1's.
            break;

        // Add 1 to reverse index, going from left to right:
        int carry = firstBit;
        // If adding 1 changed the bit to 0, we have a carry:
        while (((reverseIndex ^= carry) & carry) == 0) 
            carry >>= 1;
    }

    //combine and untangle applying phase weights
    //using successive butterflies of 2,4,...,bitwidth 
    for (int i = 1; i <= bitwidth; i++)
    {
        int halfStep = (1 << (i - 1));
        for (int j = 0; j < halfStep; j++)
        {
            int tempShift = (j << (bitwidth - i));
            w = W[tempShift];
            for (int k = j; k < length; k += halfStep * 2)
            {           
                u = complexLine[k + halfStep];
                v = complexLine[k];
                t = w * u;
                complexLine[k + halfStep] = v - t;
                complexLine[k] = v + t;
            }
        }
    }
}

//For the demonstration purposes of the project,
//we have chosen two-dimensional iterative FFT algorithm
void Solver::solve() 
{   
    /* Transform the rows */
    #pragma omp parallel for
    for (int j = 0; j < length; j++) 
    {
        FFT(complexArray[j]);
    }

    /* Transposition */
    transpose(complexArray);
    //end transposition

    /* Transform the columns */
    #pragma omp parallel for
    for (int j = 0; j < length; j++) 
    {
        FFT(complexArray[j]);
    }

    /* Transposition */
    transpose(complexArray);
    //end transposition
}
